﻿using System;
using System.Collections.ObjectModel;
using System.Globalization;

namespace BlackScholesLib
{
    internal static class NormSDist_v4
    {
        public static double NormSDist(double d)
        {
            return NormalDistribution(d, 0.0, 1.0);
        }

        /// <summary>
        /// Returns the cumulative density function evaluated at A given value.
        /// </summary>
        /// <param name="x">A position on the x-axis.</param>
        /// <param name="mean"></param>
        /// <param name="sigma"></param>
        /// <returns>The cumulative density function evaluated at <C>x</C>.</returns>
        /// <remarks>The value of the cumulative density function at A point <C>x</C> is
        /// probability that the value of A random variable having this normal density is
        /// less than or equal to <C>x</C>.
        /// </remarks>
        private static double NormalDistribution(double x, double mean, double sigma)
        {
            // This algorithm is ported from dcdflib:
            // Cody, W.D. (1993). "ALGORITHM 715: SPECFUN - A Portabel FORTRAN
            // Package of Special Function Routines and Test Drivers"
            // acm Transactions on Mathematical Software. 19, 22-32.
            int i;
            double del, xden, xnum, xsq;
            double result, ccum;
            double arg = (x - mean) / sigma;
            const double sixten = 1.60e0;
            const double sqrpi = 3.9894228040143267794e-1;
            const double thrsh = 0.66291e0;
            const double root32 = 5.656854248e0;
            const double zero = 0.0e0;
            const double min = Double.Epsilon;
            double z = arg;
            double y = Math.Abs(z);
            const double half = 0.5e0;
            const double one = 1.0e0;

            double[] a =
            {
                2.2352520354606839287e00, 1.6102823106855587881e02, 1.0676894854603709582e03,
                1.8154981253343561249e04, 6.5682337918207449113e-2
            };

            double[] b =
            {
                4.7202581904688241870e01, 9.7609855173777669322e02, 1.0260932208618978205e04,
                4.5507789335026729956e04
            };

            double[] c =
            {
                3.9894151208813466764e-1, 8.8831497943883759412e00, 9.3506656132177855979e01,
                5.9727027639480026226e02, 2.4945375852903726711e03, 6.8481904505362823326e03,
                1.1602651437647350124e04, 9.8427148383839780218e03, 1.0765576773720192317e-8
            };

            double[] d =
            {
                2.2266688044328115691e01, 2.3538790178262499861e02, 1.5193775994075548050e03,
                6.4855582982667607550e03, 1.8615571640885098091e04, 3.4900952721145977266e04,
                3.8912003286093271411e04, 1.9685429676859990727e04
            };
            double[] p =
            {
                2.1589853405795699e-1, 1.274011611602473639e-1, 2.2235277870649807e-2,
                1.421619193227893466e-3, 2.9112874951168792e-5, 2.307344176494017303e-2
            };


            double[] q =
            {
                1.28426009614491121e00, 4.68238212480865118e-1, 6.59881378689285515e-2,
                3.78239633202758244e-3, 7.29751555083966205e-5
            };
            if (y <= thrsh)
            {
                //
                // Evaluate  anorm  for  |X| <= 0.66291
                //
                xsq = zero;
                if (y > double.Epsilon) xsq = z * z;
                xnum = a[4] * xsq;
                xden = xsq;
                for (i = 0; i < 3; i++)
                {
                    xnum = (xnum + a[i]) * xsq;
                    xden = (xden + b[i]) * xsq;
                }
                result = z * (xnum + a[3]) / (xden + b[3]);
                double temp = result;
                result = half + temp;
            }

                //
            // Evaluate  anorm  for 0.66291 <= |X| <= sqrt(32)
            //
            else if (y <= root32)
            {
                xnum = c[8] * y;
                xden = y;
                for (i = 0; i < 7; i++)
                {
                    xnum = (xnum + c[i]) * y;
                    xden = (xden + d[i]) * y;
                }
                result = (xnum + c[7]) / (xden + d[7]);
                xsq = Math.Floor(y * sixten) / sixten;
                del = (y - xsq) * (y + xsq);
                result = Math.Exp(-(xsq * xsq * half)) * Math.Exp(-(del * half)) * result;
                ccum = one - result;
                if (z > zero)
                {
                    result = ccum;
                }
            }

                //
            // Evaluate  anorm  for |X| > sqrt(32)
            //
            else
            {
                xsq = one / (z * z);
                xnum = p[5] * xsq;
                xden = xsq;
                for (i = 0; i < 4; i++)
                {
                    xnum = (xnum + p[i]) * xsq;
                    xden = (xden + q[i]) * xsq;
                }
                result = xsq * (xnum + p[4]) / (xden + q[4]);
                result = (sqrpi - result) / y;
                xsq = Math.Floor(z * sixten) / sixten;
                del = (z - xsq) * (z + xsq);
                result = Math.Exp(-(xsq * xsq * half)) * Math.Exp(-(del * half)) * result;
                ccum = one - result;
                if (z > zero)
                {
                    result = ccum;
                }
            }

            if (result < min)
                result = 0.0e0;
            return result;
        }
    }
}
